As computing power of computers has improved in recent years, methods for computer simulations have developed. As a result, computer simulations have been used for various fields of application.
As numerical calculation methods for solving problems of the continuum such as in fluids and elastic bodies, the finite difference method, finite element method, or finite volume method or the like for solving the approximate solution of a differential equation on the basis of a mesh have often been used. Since numerical calculation has been utilized in fields of application such as CAE (Computer Aided Engineering), these numerical calculation methods have developed, and problems in which a fluid and a structure interact with each other have been solved.
However, in a method for using a mesh for numerical calculation, handling of cases where a moving boundary is generated is complicated, such as in a problem in which an interface such as a free surface is present and in a fluid-structure interaction problem, so program creation is often difficult.
In contrast, in a particle method (such as an MPS (Moving Particle Semi-implicit) method, SPH (Smoothed Particle Hydrodynamics) method, or the like) in which a mesh is not used for numerical calculation, special treatment is not needed in the handling of the moving boundary. Therefore, particle methods have been widely used in recent years.
FIG. 1 is a diagram explaining a problem of the particle method.
The particle method has been developed with an aim of handling easily a boundary which deforms and moves, such as a free surface or the like. Since the continuum becomes a mere particle group composed of particles 11 due to a discretization method, where the boundary surface of the continuum exists is unclear, as illustrated in FIG. 1. Therefore, there has been no unified method in the particle method for problems in which a boundary of surface tension or the like has to be visibly treated.
As for finite difference method and finite element method, there exists a calculation method for surface tension called a CSF (Continuum Surface Force) model (for example, see Non-Patent Document 1). The CSF model was applied to the particle method by X. Y. Hu, and N. A. Adams, and the calculation of surface tension by the particle method was thereby enabled (for example, see Non-Patent Document 2). In addition, there exists a method for determining surface tension as an interaction between particles (for example, see Patent Document 1), and a method for determining surface tension on the basis of deviation of neighboring particles from the position of the center of gravity (for example, see Non-Patent Documents 3, 4).
Here, an outline of the calculation method based on the CSF model applied to the particle method will be explained.
FIG. 2 is a diagram illustrating an influence radius and a continuous function obtained by overlapping kernel functions.
First, a color function <Ci> used for calculation of surface tension is defined as follows:
                              〈                      C            i                    〉                =                              ∑            j                    ⁢                                          ⁢                                                    m                j                                            ρ                j                                      ⁢                          W              ⁡                              (                                                                                                                        r                        i                                            -                                              r                        j                                                                                                  ,                  h                                )                                                                        [                  Formula          ⁢                                          ⁢          1                ]            wherein a subscript i expresses the value of i-th particle 11, and j expresses the value of j-th particle 11.r,m,ρare the position vector, mass, and density of the particle 11, respectively, and a value is given to each particle 11.Wis an overlap of kernel functions 22 used in the SPH method, and a cubic spline function 23 or the like is often used.hexpresses a radius of a spherical area in which the value of the kernel function 22 is non-zero, and is called as an influence radius 21 as illustrated in FIG. 2.
The gradient of the color function (∇Ci) is a standard method for the SPH method, and is calculated as follows.
                              ∇                      C            i                          =                              ρ            i                    ⁢                                    ∑              j                        ⁢                                                  ⁢                                                            m                  j                                (                                                                            〈                                              C                        j                                            〉                                                              ρ                      j                      2                                                        +                                                            〈                                              C                        i                                            〉                                                              ρ                      i                      2                                                                      )                            ⁢                                                ∂                                      W                    ⁡                                          (                                                                                                                                                            r                              i                                                        -                                                          r                              j                                                                                                                                ,                        h                                            )                                                                                        ∂                                      r                    i                                                                                                          [                  Formula          ⁢                                          ⁢          2                ]            
The surface tension applied to the particle i is defined as follows using the above [Formula 2].
                                          f                          s              ,              i                                =                                    ∑              j                        ⁢                                                  ⁢                                          m                i                            ⁢                                                m                  j                                ⁡                                  (                                                                                    T                        i                                            +                                              T                        j                                                                                                            ρ                        i                                            ⁢                                              ρ                        j                                                                              )                                            ⁢                                                                    r                    i                                    -                                      r                    j                                                                                                                                r                      i                                        -                                          r                      j                                                                                                    ⁢                                                W                  ′                                ⁡                                  (                                                                                                                                    r                          i                                                -                                                  r                          j                                                                                                            ,                    h                                    )                                                                    ⁢                                  ⁢        wherein        ⁢                                  ⁢                  T          i                =                              Γ            g                    ⁢                      1                                                        ∇                                  C                  i                                                                            ⁢                                    (                                                                    1                    d                                    ⁢                  I                  ⁢                                                                                                          ∇                                                  C                          i                                                                                                            2                                                  -                                                      ∇                                          C                      i                                                        ⊗                                      ∇                                          C                      i                                                                                  )                        .                                              [                  Formula          ⁢                                          ⁢          3                ]            Γg is a surface tension coefficient,dis the dimension(s), andIis an identity matrix.
The calculation method based on the CSF model is a calculation method for performing time development using the surface tension obtained by the above [Formula 3].
FIGS. 3-7 are diagrams illustrating the temporal changes of the result of using the conventional calculation method based on the CSF model.
As illustrated in FIGS. 3-7, in the calculation method based on the CSF model, a fluid parcel falls due to gravity receiving the influence of surface tension, and collides against a plane surface. The calculation method disclosed in Non-patent Document 2 is used for the calculation of the behaviors illustrated in FIGS. 3-7. In FIGS. 3-7, the vertical direction is set to be the y-axis so that y=0 expresses a wall surface, and the fluid 31 cannot enter the position where y<0.
As illustrated in FIG. 3, the initial shape of the fluid 31 is a square of 0.05[m]×0.05[m]. In addition,
Γg=10
[N/m2],
ρi=1000
[kg/m3],
h=0.0046875
[m],
c=109.8
[m/s],
dt=1.423×10−6 
[s], viscosity
μ
is 0.001[m2/s], gravity acceleration
g=9.8
[m2/s], and a cubic spline function is used as a kernel function. The particles in the initial state are positioned in a uniform grid pattern (33 particles per side, for a total of 1,089 particles), and the velocity of all particles is zero.
As for the calculation results, as illustrated in FIG. 4, the distribution of the fluid 31 is in a smooth circular shape at time 0.04[s], and as illustrated in FIG. 5, the fluid 31 collides against the wall at time 0.14[s] and greatly deforms.    [Patent Document 1] Japanese Laid-open Patent Publication No. 2008-111675    [Non-Patent Document 1] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti, “Modelling Merging and Fragmentation in Multiphase Flows with SuRFER”, Journal of Computational Physics, 113, 134, 147, (1994)    [Non-Patent Document 2] X. Y. Hu, N. A. Adams, “A multi-phase SPH method for macroscopic and mesoscopic flows”, Journal of Computational Physics, 213, pp. 844-861 (2006)    [Non-Patent Document 3] T. Hongo, M. Shigeta, S. Izawa, and Y. Fukunishi, “Sanjigen Hiassuku SPH Hou ni okeru Kiekikaimen ni Sayousuru Hyoumenchouryokumoderu (Surface Tension Model Acting on Gas-Liquid Interface in Three-Dimensional Incompressible SPH Method)” Dai 23 kai Suuchi Ryuutairikigaku shinpojiumu (The 23rd Numerical Fluid Dynamics Symposium), A8-5 (2009)    [Non-Patent Document 4] M. Agawa, M. Shigeta, S. Izawa, and Y. Fukunishi, “Shamenjyo wo Ryuugesuru Ekitai no Hiasshuku SPH Simureshon (Incompressible SPH Simulation of Fluid Flowing Down on Slope)”, Dai 23 kai Suuchi Ryuutairikigaku shinpojiumu (The 23rd Numerical Fluid Dynamics Symposium), A9-4 (2009)